#ifndef __FUNCTION__
#define __FUNCTION__

#include <iostream>
#include <cmath>
const double eps = 1e-12;

//根据课上所讲，设计了一个函数类，利用了纯虚函数及相关的函数必备，如导数，求值。
class Function {
public:
    virtual double operator()(const double& x) const = 0;
    virtual double get_diff(const double& x) const { 
		return ((*this)(x + eps) - (*this)(x - eps)) / (2 * eps); 
	}
    virtual double get_secdiff(const double& x) const { 
		return ((*this)(x + 2 * eps) + (*this)(x - 2 * eps) - 2 * (*this)(x)) / (4 * eps * eps); 
	};
};


// 多项式为函数的子类，我们本次大作业也主要是进行相应的多项式生成
class Polynomial : public Function{

public:
    // 基本的构造函数
    Polynomial() : coef({ 0 }), n(0){};
    // 常数多项式
    Polynomial(const double c0) : coef({ c0 }), n(0){};
    // 一次多项式
    Polynomial(const double c0, const double c1) : coef({ c0, c1 }), n(1){};
    // 两次函数
    Polynomial(const std::vector<double>& _coef) : coef(_coef), n(_coef.size() - 1){};
    //析构函数要释放向量空间coef
    ~Polynomial() {
     coef.clear();
    }

    //求值，对每个次数分别求
    double operator()(const double& x) const {
        double res = 0, x_degree = 1;
        for (int i = 0; i <= n; i++) {
            res += coef[i] * x_degree;
            x_degree *= x;
        }
        return res;
    };

    //求导
	double get_diff(const double& x) const {
		double res = 0, x_degree = 1;
		for (int i = 1; i <= n; i++) {
			res += coef[i] * x_degree * i;
			x_degree *= x;
		}
		return res;
	}
    //二阶导数，本次作业没用到
	double get_secdiff(const double& x) const {
		double res = 0, x_degree = 1;
		for (int i = 2; i <= n; i++) {
			res += coef[i] * x_degree * i * (i - 1);
			x_degree *= x;
		}
		return res;
	}

    //下面分别定义了多项式的加减乘
    Polynomial operator+(const Polynomial& rhs) {
        Polynomial res = *this;
        for (int i = 0; i <= rhs.n; i++) {
            if (i <= res.n) {
                res.coef[i] += rhs.coef[i];
            } else {
                res.coef.push_back(rhs.coef[i]);
                res.n++;
            } 
        }
        return res;
    };
    Polynomial operator-(const Polynomial& rhs) {
        Polynomial res = *this;
        for (int i = 0; i <= rhs.n; i++) {
            if (i <= res.n) {
                res.coef[i] -= rhs.coef[i];
            } else {
                res.coef.push_back(-rhs.coef[i]);
                res.n++;
            } 
        }
        while (res.coef[res.n] == 0 && res.n > 0) {
            res.coef.pop_back();//因为多项式私有类有其次数，所以减法要考虑其降次
            res.n--;
        }
        return res;
    };
    Polynomial operator*(const Polynomial& rhs) {
        Polynomial res;
        res.n = n + rhs.n;
        res.coef.resize(res.n + 1);
        for (int i = 0; i <= n; i++)
            for (int j = 0; j <= rhs.n; j++) {
                res.coef[i + j] += coef[i] * rhs.coef[j];
            }

        return res;
    };
    Polynomial diff() const {
        Polynomial res;
        res.n = n - 1;
        res.coef.resize(res.n + 1);
        for (int i = 0; i <= res.n; i++) {
            res.coef[i] = coef[i + 1] * (i + 1);
        }
        return res;
    };

    // 输出函数，将多项式输出,考虑到要用matlab画图，将形式也控制一下
    friend std::ostream& operator<<(std::ostream& os, const Polynomial& poly) {
        std::cout << poly.coef[0] << std::endl;
        if (poly.n == 0) {
            return os;
        } else {
            for (int i = 0; i <= poly.n; i++) {
                if (poly.coef[i] == 0) continue;
                if (poly.coef[i] > 0) {
                    std::cout << " + " << poly.coef[i] << " * x.^" << i;
                } else {
                    std::cout << poly.coef[i] << " * x.^" << i;
                }
            }
            std::cout << std::endl;
            return os;
        }
    };
    private:
    int n;      // 多项式的次数
    std::vector<double> coef;       // 多项式次数
};


// 因为B样条是确定了节点就相应固定，所以这里作为函数子类声明。
class BsplineBase : public Function{
private:
    int N;//度数
    std::vector<double> knots;// 插值的节点

public: 
    BsplineBase(int _N, std::vector<double> _knots) : N(_N), knots(_knots) {
        if (knots.size() != N + 2) {
	        std::cerr << "度数不匹配" << std::endl;
            exit(-1);
        }
    }
    BsplineBase(std::vector<double> _knots) : N(_knots.size() - 2), knots(_knots) {
        if (knots.size() < 3) {
            std::cerr << "输入数据错误" << std::endl;
            exit(-1);
        }
    }
    //析构函数
    ~BsplineBase() {
        knots.clear();
    }
    //这里就是按照课本定义，从一次开始生成
    double operator()(const double &x) const {
        if(N == 1) {
		    if(x > knots[0] && x < knots[1])
		        return (x - knots[0]) / (knots[1] - knots[0]);
		    else if(x >= knots[1] && x < knots[2])
		        return (knots[2] - x) / (knots[2] - knots[1]);
		    else
		        return 0;
	    } else {
            //生成完一次样条后就可以进行迭代
            
	        std::vector<double> knots1(knots.begin(), knots.end() - 1);
	        std::vector<double> knots2(knots.begin() + 1, knots.end());
	        BsplineBase b1(N - 1, knots1);
	        BsplineBase b2(N - 1, knots2);
	        double c1 = (x - knots[0]) / (knots[N] - knots[0]);
	        double c2 = (knots[N+1] - x) / (knots[N+1] - knots[1]);
	        return c1 * b1(x) + c2 * b2(x);
        }
    }
};


#endif
